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OF ATMOSPHERIC GUST RESPONSE PROBLEM 
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SUMMAEX 


Equations of motion are derived for the dynamic response of an air- 
craft to random atmospheric gust loads. These equations include the 
degrees of freedom of plunging, pitching, rolling, and an arbitrary num- 
ber of elastic normal modes. Solutions of these equations are expressed 
in terms of a number of so-called primitive solutions obtainable by intro- 
ducing the Dirac delta function. The solutions for center -of- gravity 
acceleration response and wing-root bending-moment response depend upon 
certain autocorrelations and cross correlations which enter the analysis. 
Results for simplified cases show that unsteady aerodynamic theory is 
not important for increasingly large values of the turbulence scale com- 
pared with values of the wing chord. However, the pitching degree of 
freedom exhibits an important effect as the turbulence scale increases. 

The results are also compared with the results of the usual sharp-edged- 
gust formula. 


INTRODUCTION 


The subject of aircraft gust loads has been studied by theoreticians 
anri experimentalists for many years. It is one of the principal consid- 
erations in the design of both military and commercial aircraft; that is, 
airworthiness requirements for all types of aircraft include gust-load 
criteria. These criteria in all cases are based on the sharp- edged-gust 
formula developed by Rhode in 1931 (ref. 1) and in most cases take into 
account the effects of gust alleviation due to rigid-body motion. 

The engineer's understanding of the gust-load problem has steadily 
increased over the years through basic research investigation and opera- 
tional experience. Foremost among these investigations is the work carried 
out by the National Advisory Committee for Aeronautics. A great body of 
exper ime ntal data Is now available to engineers and aircraft designers 
through these research efforts. Important statistical information relating 
to peak gust loads experienced by aircraft of various types subject to 
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atmospheric gusts in various types of atmospheric and topographical 
environment has been collected by VGH recorders (ref. 2). Such informa- 
tion may be used to determine the probability density distribution of 
peak gust loads as a f motion of airspeed, altitude, and several important 
aircraft parameters. 

More recently, important information has become available on the 
gust structure in the atmosphere through the application of generalized 
harmonic analysis (refs. 3 and 4). This type of analysis has been used 
to predict gust loads. In addition to these analyses, which are in the 
nature of applied research, more basic theoretical investigations which 
bear directly on the gust-loads problem may be cited. For example, 
recent papers by Liepmann (refs. 3 and 6) and Fung (ref. 7) are among 
the important contributions to the present state of knowledge. 

A fundamental aspect of a VI of these investigations is the necessity 
for giving the gust structure a mathematical form convenient for inclusion 
in a rational analysis of gust loads. Unfortunately, this is still a weak 
area which can be strengthened only by continued experimental and theoret- 
ical investigation. In this connection effective use of turbulence theo- 
ries may prove to be very fruitful. Certainly, the problem of atmospheric 
turbulence bears too close a resemblance to many of the turbulence prob- 
lems in aeronautics already studied to be ignored. It is well known, of 
course, that many of these latter “classical" problems are far from being 
solved. Indeed, the subject of turbulence in general has stimulated some 
of the keenest intellects in historical times, yet many areaB of the tur- 
bulence problem remain unconquered. 

In the present report the concepts derived by Liepmann and Fung are 
extended to include the complications introduced by a multi-degree- of - 
freedom system, the airplane. In some respects the analysis parallels 
the work of Diederich (ref. 8) . Related experimental investigations are 
reported in references 9 and 10. 

This work was conducted at the Massachusetts Institute of Technology 
under the sponsorship and with the financial assistance of the NACA. 


SYMBOLS 


A r (x,y) 

b 


b c 

C(k) 


mode shape of rth normal mode, positive up 
wing semiBpan 
wing semichord 
Theodor sen function 
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cl 

cja(y) 

c (y) 

Ee 

E± 

[•] 

[ p jj] 

f jjW 

J 

k 

L 

L(x,y,t) 

l 

V 


M 

m a 

% 

Mr 

Mp*(cn) 


lift coefficient 

section lift-curve slope of wing at station y 
chord of wing at station y 
external work of applied forces 
internal strain energy 
generalized stiffness matrix 
forcing-function matrix 

diagonal element of [Fj 
system function 

pitch radius of gyration referred to wing semichord 
reduced frequency, aib c /U 
integral scale of turbulence 
lift on the wing at point (x,y) 

distance traveled in arbitrary time interval, Ut 

tail length measured from airplane center of gravity 

distance between wing midchord and tail midchord in plane 
of symmetry 

total number of modes 

mass of airplane 

bending moment 


generalized mass. 



dm A r (x,y) 


transfer function 
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mass matrix 

N number of generalized coordinates 

P velocity potential 

p tt d 2 P/St 2 

p pressure 

Ap(x,y^t) pressure difference at point (x,y), positive when it pro- 
duces lift 

q generalized coordinate 

q dynamic pressure 

R ratio of horizontal tail area to wing area 

r, s normal modes (subscript or superscript) 

S total wing area 

S' projection of airplane plan form in xy-plane 

T kinetic energy of airplane 

t time 

U mean forward speed of airplane 


V=» z + x0 + AfCx^yOij 


v 

W 

w 

v* 

x,y 


croBsflow 

stochastic function with zero mean value 
gust velocity 

root-mean-square intensity of turbulence 

coordinate system, x positive, forward, y positive along 
right wing 
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xo>y 

x,y 

#w 

x^y 

z 


a 



N 


8 

8 e 

S 

T] 

©r(x,y) 

0 

K 

Aj(y) 

^r>^s 

M- 

V 

£r> £s 

P 

a 


coordinates of midchord 
coordinates of three-quarter chord 
coordinates of quarter chord 

vertical displacement of center of gravity of rigid air- 
plane relative to average position, positive up 

angle of attack 

equivalent gust angle of attack 

aerodynamic matrices 
sweep angle 
unspecified matrix 

unit impulse 
control angle 

integration variable (streamwise) 
integration variable (spanwi.se) 

mode shape for twisted surface which is rigid chordwise 

angular displacement about center of gravity of rigid air- 
plane, positive nose up 

airplane relative density 

spanwise weighting function 

bending -moment coefficients for r,s normal modes 

generalized dependent variable 
integration variable (streamwise) 

displacement of r,s normal modes, positive up 
air mass density 
standard deviation 
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integration variable (time) 
power spectral density- 

bank angle of rigid airplane , positive right wing up 
correlation function 
correlation interval (distance) 

Kussner function 
angular frequency 

natural (angular) frequency of rth normal mode 

derivatives with respect to time 

« 

derivatives with respect to distance (measured along 
x-axis) 

average 

DERIVATION OF DYNAMICAL EQUATIONS OF MOTION 


The problem of aircraft response to atmospheric turbulence (gust 
response) will be considered in a rather general manner in order to bring 
out some of the important features of the problem. The analysis is set 
up on the basis of familiar dynamical principles and is carried through 
as far as possible without introducing statistical concepts. However, 
a certain point is reached in the analysis, when solutions to the equa- 
tions of motion are desired, where the introduction of statistical con- 
cepts becomes mandatory if the analysis is to continue. These concepts 
are then brought into the analysis and the final solution iB obtained. 

The analysis Is carried out this way in order to show that a solution to 
the problem depends upon the proper use of familiar dynamical analysis 
principles as well as upon statistical analytical tools . 

The problem may be formulated by requiring that the "response" of an 
aircraft flying with velocity U through turbulent air is to be deter- 
mined. The particular response of interest may be the acceleration, veloc- 
ity, displacement, bending moment, shear, stress, and so forth at any 
arbitrary point of the aircraft structure. 

The analysis will include three rigid-body degrees of freedom, heaving, 
pitching, and rolling (z, 9, and 0, respectively), and an unspecified 

number of elastic degrees of freedom (| r where r = 1, 2, 3, . . .) 
defined by the principal vibration modes of the aircraft; that is, the 
analy sis Is constructed for normal modes of the entire aircraft defined 
by A r (x,y) (fig. 1). 
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Use is made of Lagrange's equation to establish the equations of 
motion for the system, for small perturbations from eq uil ibrium, namely 
(ref. 11) : 


d/^T \ t dEj _ 5E e 
dtldqJ dq dq 


The kinetic energy is 


= i T v 2 dm 

2 J s « 


Thus, 


■where 


_d/ST 
d? 


r V^am 

dt Jgi 3S[ 


V=z + x0+y#+^ A r (x,y)| 3 


The internal strain energy is (ref. 12) 


where 




3q 


Mrqr 2 lr 


(9. 4 £r) 

Co. = l r ) 
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The potential energy of the external forces is 


Ee = 



z + 0x + + 


£W>! 


f!p(x,y,t) dy dx 


The pressure difference acting on the aircraft hp(x,y,t) contains two 
distinct contributions. The first part is due to disturbances of the 
field caused directly by motions of the aircraft. This part may be cor- 
rectly accounted for by solution of the wave equation with appropriate 
boundary conditions (ref. 12) 


- p tt 


The other part is due to disturbances originating elsewhere in the fluid 
(usually at great distance from the aircraft) and, for incompressible 
flow, can be described by the Navier-Stokes equations (ref. 13, p. 215) 




.2 
8 Uj 


There may appear to be a contradiction in the above discussion, 
namely, that motion effects are determined from perfect-fluid considera- 
tions through the concept of a velocity potential, while the disturbance 
effects may be due to turbulence, which, in general, cannot be described 
by a velocity potential. However, the scale of the turbulence which is 
important for the aircraft gust problem is large enough so that viscous 
effects may be neglected; that is, viscosity is associated with the very 
small eddies in the fluid. With viscosity thus eliminated, there is no 
contradiction. 

Solution of the problem as outlined for the case of a turbulent flow 
superposed on uniform potential flow is an extremely difficult undertaking 
even in incompressible flow. In order to gain further insight into the 
problem, therefore, it seems highly desirable to make some simplifying 
assumptions. 

Of all the parameters which enter the Navier-Stokes equation, namely, 
pressure, density, velocity, and viscosity, it may be assumed that, so 
far as the present application is concerned, velocity plays the dominant 
role. Furthermore, since interest is centered upon aircraft response in 
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the three rigid-body modes of heaving, pitching, and rolling and in a 
number of elastic modes, all components of fluid velocity except that 
which is par all el to the initial z-axis may be safely discarded. Thus, 
the significant part of the problem reduces to that of finding a satisfac- 
tory description of one component of velocity, namely, w. For small per- 
turbations one can consider an equivalent angle of attack w/U = o^. The 

reasoning here can be made general in the sense that it can be applied to 
incompressible and compressible flow alike. It should be mentioned that 
the elimination of the effect of the static-pressure influence in the 
disturbance field is not completely straightforward, but it simplifies 
the problem considerably. 

The influence of the motion of the aircraft, as mentioned previously, 
can be handled as a straightforward problem in unsteady potential flow; 
this is very tedious, however, especially if al 1 boundary conditions are 
correctly included in the analysis. It is usual to make some simplifying 
assumptions here also. 

If it is decided that structural response is of primary interest, 
for example, bending moment or shear, it is convenient to consider sim- 
plifying the aerodynamics picture as much as is consistent with the desired 
accuracy. This is particularly desirable in a statistical analysis, which 
can become unwieldy for even the simplest of problems. 

To show how the statistical picture can get complicated, consider 
the factors that would, be introduced by including an additional component 
of velocity in the analysis, say the crossflow v. This means that instead 
of "one-component" turbulence as given by the perturbation w, another 
vector quantity v has been added and accordingly "two -component" turbu- 
lence must be considered. This means that additional, statistical quanti- 
ties are brought into the analysis. Important statistical parameters such 
as probability functions, correlation functions, and spectra would, have 
to reflect this fact. In the general case of two-component turbulence, 
the correlation tensor Rji(s) would contain four terms. This may be 

compared with a single quantity Rn(s) in the relatively simple case 

of one-component turbulence. The probability functions for two-component 
turbulence are correspondingly more difficult to determine, and the spectra 
become extremely difficult to calculate in the general case. 

These reflections add considerable weight to the argument for sim- 
plification of the analysis. Indeed, it is a matter of practical interest 
that such simplifications are a necessary adjunct to keeping the mathematics 
within reasonable proportions. 


•'■Note that a small effect may be expected because of crossflow v 
arising on the vertical tail and a somewhat larger effect of crossflow 
for swept wings may be expected. However, these effects of crossflow 
are neglected. 
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The formulation of the equations from this point requires a specifica- 
tion with regard to the aerodynamic theory ; namely, will the unperturbed 
(steady) flow conditions be incompressible or compressible and, if the 
latter, what Mach number characterizes the flow? Also, how does one 
include the aspect-ratio effects? There are a number of theories which 
afford an approximate description of 2p(x,y,t) . The theory which is 
vised will depend largely upon the free-flight Mach number and the plan 
form of the wing. In order that the present analysis may have a definite 
framework within which to proceed, the following specifications are made 
regarding the aerodynamic theory: 

(1) Aerodynamic strip theory is used for both motion effects and 

gust effects 

(2) Quasi-steady aerodynamic theory for incompressible flow is used 

to predict all motion forces by resolving the angle of attack 
at 3/4 chord 

(3) Unsteady aerodynamic theory for incompressible flow is used to 

predict gust forces (later in the analysis, this will be 
modified) 

(4) For simplicity, only the forces on the wing and horizontal tail 

are considered in the analysis 

According to the first statement, the analysis must be restricted 
to aircraft of "large" aspect ratio, but it may also include swept as 
well as straight wings. The second statement restricts the analysis to 
flight configurations where all free-body motions ( including the vibratory 
modes) occur at low values of reduced frequency, namely, k < 0 . 03 ; and , 
of course. Incompressible flow requires that the Mach n umb er be low. 


On the basis of these assumptions Eg may be written as follows : 




3 v/'b £* + + + ^ A r(x,y)| r 

J t \z + ex’ + 0y +^A r (x',y)i 2 


L(y>x,t) dy + 

dy 


The lift forces on the wing L(y,x,t) dy and on the tail L*(y,x',t) dy 
are assumed to be concentrated at the local quarter-chord points of t ho se 
surfaces. 
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There are two parts to L(y,x,t) and L'(y,x',t), corresponding to 
the forces due to motion and the forces due to turbulence. The forces 
due to motion may be calculated by resolving the angle of attack at the 
three-quarter-chord point as follows (ref. 12) : 2 


L(y^t) m = qc(y)C z 


a 


r -1 


e + 


^9r(x,y)£ r - 


z + 0x + #y + \ ^(5^)1 


L 1 


u 


l* , (yj‘3E , ,t) m — g , c , (y)c l ' 

w 


z + 0x‘ + #y + YA r (x , ,y)£ r 
9 + /_ j 0r ^ 5,y ^ r 


U* 


The dynamic pressure is given the symbol q to distinguish it from the 
generalized coordinate q. It is important to keep account of the coordi- 
nates involved in the above expressions; unprimed quantities refer to 
wing coordinates while primed quantities refer to tail coordinates, x,y 
is the local coordinate of the wing quarter chord, and x,y is the local 
coordinate of the wing three-quarter chord. Note that a new quantity is 
introduced in addition to the new coordinates, that is, the qu ant ity 0 r 
which gives the elastic twist of the wing (or tail) at the three- quarter- 
chord point. 

The other part of the aerodynamic force, that part due to turbulence, 
may be calculated from the relations (ref. 12) 


L(y,X,t) g => qc(y)C 2 a F *(t - T,y)a i (T,y) dT 

U -CO 

r” * 

L*(y,2',t) g = q' c ' (y)C j 1 / *'(t - r^oj^y) dT 

^ —00 


%ote that effect of the wing wake vorticity on the horizontal tail 
is here neglected for simplicity since it is debatable whether much is 
to be gained by including this effect in an analysis of this complexity. 
Probably the only time it can be fully justified is when tail loads or 
tail stresses are of primary interest. 
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In these relations the specifications of the kernel function t are 
purposely left in general terms in order to avoid confusion. However, it 
may be said that this function must take into account sweep, taper, and 
the finite time required for the gust pattern to pass from wing to tail. 

In this connection it should be brought out that Taylor's hypothesis is 
assumed valid for this analysis; that is, the gust field is regarded as 
a "frozen" pattern which does not alter with time locally, but rather 
the time dependence is caused by the uniform motion (at velocity U) of 
the aircraft. This may be expressed in succinct terms by stating merely 
that S/St =0 in the fluid. 

This completes the formulation of the problem sufficiently to allow 
derivation of the dynamical equations of motion for the aircraft. Direct 
substitution of the above results into Lagrange's equation yields the fol- 
lowing matrix equations: 

H {«} + H 10 + H W + H W - Ph] W 

The various square matrices are given in the appendix. The generalized 

coordinate matrices jqj contain the generalized coordinates z, 0, 

0 , and § r (r = 1 , 2 , 3, . . . ) . It is interesting to note in passing 

that the following aeroelastic problems may be dropped out of the above 
equation as special cases: Divergence, 

[a] {q} + [e] {q} = 0 (l) 

stability, 

[m]{*q‘] +[S]{q} + [a] {q} + [«] { <1 } =° ( 2 ) 

control response, 

[m]{qj +[a]{qj + [a] jq j +[e]jqj = [AlJ { l] B e (t) (3) 

and flutter. 


-u£[m'] | q | + icoC<k)[a] j^qj + C(k) [a] jqj + [e] jqj = 0 (4) 
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(The matrix [m'J includes apparent mass effects.) A slight revision 
is necessary to obtain the flutter case; namely, the Daeodorsen function 
C(k) must he included in the motion forces. For simplicity, C(k) may 
he assumed to he independent of the span coordinate. 

It is interesting to note that equation ( 3 ) provides the basis for 
studying a gust alleviator system. If QFj jj (l) is added to the right 

side of equation (3), an investigation can he made of the effectiveness 
of various systems whose behavior is governed by a relation of the form 


5 e (t) = / hg(t - t)w(t) dT 

'J -00 


where hg(t) is the system function. An interesting problem may be posed 
by requiring that h s (t) be selected for minimu m root-mean-square response 

of, say, pitching displacement. This is the so-called Wiener-Hopf type 
of problem. 


SOLUTION OF DYNAMICAL EQUATIONS OF MOTION AND 
DETERMINATION OF ROOT-MEAN-SQUAEE RESPONSE 


It was shown in the previous discussion that the equations of motion 
for an aircraft subject to turbulent velocity fluctuations in the xy-plane 
and perpendicular to the xy-plane are formulated as follows: 


[m] {’ij + [8] {4} + [a] {q} + [e] {q} - 


where the various matrices are given in the appendix. It is shown that 
the diagonal matrix jj is the forcing-function matrix and is of the 
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Note that it is a function of the parameter £ = Ufc. Various methods of 
solving this set of equations are now considered. 

It is well to point out that in the present problem fn( £) , 
fp?( O ♦ • • • are Stochastic functions of £, and solutions of the equa- 

tions must be approached in light of this fact. This increases the com- 
plexity of the problem considerably. 

For purposes of comparison, some more-or-less standard methods for 
solving linear equations are listed as follows: (l) By Laplace transforms, 

(2) by Fourier transforms, ( 3 ) by finite differences, and (4) by superposi- 
tion of indicial solutions. 

It is evident that direct application of Laplace transforms must fail 
because of the nature of Fj that is, the matrix is made up of random 

functions which are considered to have a beginning at -00 in time and 
the ini tial conditions are indeterminant. Direct application of the 
Fourier transform, even in the bilateral case, must also fail because, 
in general, the Fourier transform for the functions entering the right-hand 
side of the equation does not exist; that is, the integral may be diver- 
gent. Aside from these difficulties there is no generally consistent way 
of specifying the functions fjj(£) because they depend upon the random 
variable a i (x,y) . Even if it were possible to specify these functions 

in detail, direct application of the two transform methods would yield 
nothing because of the aforementioned difficulties. Solutions by fi n ite - 
difference techniques could, in principle, be carried out but, in practice, 
would be a formidable undertaking. This leaves the last method listed 
above as a logical choice, if not, indeed, the only choice. It will be 
shown presently that the method of superposition can be applied to deter- 
mine a limited amount of information about the problem. 

The method of superposition consists in building up a complete solu- 
tion from certain special types of solutions _which are relatively easy 
to obtain. There are two common methods available, either of which may 
be applied here, namely, superposition of simple -harmonic solutions through 
the use of the Fourier integral or superposition of unit- impulse-response 
solutions through the convolution integral. The latter seems to be the 
better method because it brings in the statistical parameters associated 
with the problem in a more straightforward manner. However, it is inter- 
esting to point out that the solutions obtained by the unit-impulse method 
are related to the solutions obtained by the simple -harmonic method through 
the Fourier transform, just as the unit-impulse response is related to the 
simple-harmonic response. 
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In considering a solution to the matrix equation by the method of 
superposition it is convenient to introduce the Dirac delta function to 
define a "primitive" solution, which corresponds to replacing the right- 
hand side of the equation by the following matrix: 




The unit impulse 8(£) is the forcing function in this instance 
and its position in the matrix is given by the subscript j . A solution 
for this case is called q* and q may be any one of the several gen- 
eralized coordinates z, 8, 0, or £ r . Clearly, then, the solution of 
interest may be written (ref. 12, p. 2 lj-5) 



f ijd - E)fjj(E) « 

_oo 


where N is the total number of generalized coordinates and l is in 
the same units as £. By definition qj(x) = 0 for x < 0. The cal- 
culation concerning the quantities qj may be carried out by any con- 
venient method, including Laplace transforms and finite differences. 
However, these calculations are especially tedious if N is much larger 
than 3 or For cases where a large number of generalized coordinates 
is used in the analysis, the use of high-speed computing machines is 
strongly recommended. 

The problem is far from being solved, however, with the writing down 
of the above equation. At this point in the analysis it is realized that 
one must apply statistical methods to learn additional information about 
q(l), which is a stochastic function simply because the quantities 

are stochastic functions . Another question must also be settled at this 
point in the analysis j namely, what is the quantity of primary interest? 
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If the quantity of interest is one of the generalized coor dinat es 
or the nth derivative of one of these coordinates, then most of the 
information is available from the equation as it stands. On the other 
hand, if a quantity such as bending moment or shear at some point in the 
wing is of interest, then additional considerations must be given to the 
solution in its present form. This point may be made clear if one con- 
siders the following two typical cases, each. one or both of which may be 
of practical interest: (l) Study of aircraft center-of-gravity accelera- 

tion, and (2) study of wing-root bending moment. 

In the first case the solution to be studied takes the form 


N 

ZjU - £) f jj(£) 

1 

while in the second case the solution to be studied takes the following 
form: 


M B (i) 


r b 

= - y dm 

^ 0 

•• ^ r— ] 

v4> + x0 + z + \ A(x,y)i r 

It 1 


r 


+ 



dy L(x,y,t) m + 


L 


L(x,y,t) g jy 

•• •• 

Notice that each of the quantities (, 0, and so forth as well as 

L(y,x,t) m are complicated functions of the corresponding primitive 

• • •• 

solutions 0j, z*j, and so forth and the forcing functions fjj 


It is apparent that the calculation of the bending moment in the 
manner indicated above is quite involved. Ah' alternative means is avail- 
able which is somewhat simpler in mathematical form; that is, the bending 
moment can be obtained by linear superposition of the normal modes as 
follows (ref. 12) : 


M B (Z) 


M 

= V ^rlr(^) 
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where the quantity A r gives the root bending-moment coefficient due 

unit deformation of the rth normal mode. It can be computed from the 
inertia loading on the wing corresponding to the rth normal vibration 
mode. The quantity M is the total number of modes that are used in 
the analysis; for this case, M = N - 3 - The quantities l r can be 

obtained by superposition of the primitive solutions discussed 

previous ly . Thus , 


M K 

M B (Z)=£ Sr jd - 5)fjj(S) 45 

x = i i=i 


to 


Consideration is given now to the determination of certain statistical 
quantities associated with the solutions of the dynamical equations of 
motion. Interest is here centered on the study of center -of -gravity 
acceleration and wing-root bending moment, respectively, as given by the 
relations 


N 




*j=*l 


N M 
j=l r=l 


Both of the physical quantities on the left of the above equations 
are, of course, stochastic functions and the methods of statistical anal- 
ysis must accordingly be applied to obtain additionally desired informa- 
tion. There are several questions one might pose at this point in regard 
to the statistical properties of K and Mg. Some of these are worth 

noting as follows: (l) What is the probability density distribution? 

( 2 ) What is the correlation function? (3) What is the power density 
spectrum? 

The probability density distribution provides information on the 
amplitude of the stochastic function. It can be shown that both z 
and Mb are Gauss ianly distributed provided that f j j is also Gausslanly 

distributed. This follows as a consequence of the central limit theorem 
and the properties of linear systems. It may be assumed that all of the 
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quantities concerned, here are Gaussian for convenience, but it is not 
a necessary assumption in what is to follow. There are, however, strong 
motivations for making this assumption; that is, (l) the Gaussian curve 
has been exhaustively studied, and accordingly much is known about it, 

(2) it may be characterized by a single parameter, namely, the mean sq uar e 
or standard deviation, ( 3 ) the characterization by the mean- sq uar e param- 
eter provides a link between amplitude spectra and power spectra, and 
(4) many physical quantities follow the Gaussian distribution. 

The correlation function and power density spectrum are Fourier 
transform pairs. The correlation function provides information on the 
coherence of the stochastic function at adjacent intervals in time and/or 
space, while the power density spectrum gives information as to the dis- 
tribution of power with frequency or wavelength. The use of the word 
power is rather generalized in that it need not refer to actual physical 
power. Perhaps a better terminology would be the mean-square density 
spectrum since the area under the power spectral density curve is equal 
to the mean square of the stochastic function. However, the term "power 
spectrum" seems to be rather firmly entrenched in the literature. 

If the probability density of a stochastic function W is called 
p(W), the autocorrelation function of W Is called cp w (x), and the 

power density spectrum of W is called <& w ( cd) , then one may write 
(ref. 14) 


W 


cj2 = P w 2 p(w) dW = r & w (cd) ^ = qw 

—00 u —00 


( 0 ) 


which means that the mean-square value is obtainable from all three of the 
statistical quantities discussed above. In the following discussion, 
it is assumed that W = 0; that is, the mean value of n~l 1 stochastic 
functions is assumed to be zero. Also, the analysis is directed toward 
the determination of autocorrelation functions for (l) center -of -gravity 
acceleration and (2) wing-root bending moment. From these results it is 
possible to compute the root-mean-square response. 


Autocorrelation of Center-of -Gravity Acceleration 

The autocorrelation of the center-of-gravity acceleration may be 
formed from the solutions to the dynamical equations of motion as follows: 
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2J(l)z(l + x) = 


N N 


II L dv + x - v) x 


t=l j=l 


f d E 2 3 (£)^ j (2 - E) 

O -CO 


-1 av 


where the bar over the symbols on the left-hand side and the brackets 

on the right-hand side indicate that the average is being taken 
av 

with respect to the ensemble at a particular value of l. 

With an appropriate change of independent variable and interchanging 
the order of integration and averaging operations , the autocorrelation 
of Jz may be written 

N N 

z(l)z(l + x) = ^ 

1 j 

In order to reduce the autocorrelation of z to this form it is also 
necessary to assume that the forcing functions f± and f j are stationary 

random functions; that is, the average product f pf j must be independent 

of parallel translations of these functions along the abscissa axis. 

Another way of saying the same thing in regard to z is to say that the 
left-hand side of the equation is independent of l for a stationary 
process. (These remarks do not mean, however, that the functions are 
ergodic . ) 

The double integral on the right may be reduced to a single integral 
by defining the new quantity 


jT &i(£)zj(£ - vJf^Cx - v) d£ dv 


C ^0 pOO 

q>ij (-v) = / - v) d£ 

u —00 


which is the cross correlation between the ith primitive solution and 
the jth primitive solution corresponding to In general, there will 

be such cross-correlation functions needed in the analysis. These 
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quantities depend, on the primitive solutions obtained by successively 
substituting -unit-impulse functions on the right-hand side of the dynami- 
cal equations of motion. These cross-correlation functions are referred 
to as the primitive cross -correlation functions corresponding to z. 

After substitution, the autocorrelation for z reads 


N N 


2z(x) " I I cpij^C-v) fff j(x - v) dv 


There remains the problem of computing fjfj, but before this is 

considered the autocorrelation of wing-root bending moment will be 
discussed. 


Autocorrelation of Wing-Root Bending Moment 

The autocorrelation of wing-root bending moment may be formulated 
as follows: 


N N M M 


M B (2)M B (2 + x) = Y, Z X X ^ IT ” v ) f i f j( 5 " v ) dv d C 


i j r s 


In a manner entirely analogous to the formulation of 2 autocorrela- 
tion, a new cross-correlation function may be defined as 


9ij S (-v) = J 6ri(£)6 B jC5 “ v) 

In general, there will be (NM) 2 such cross-correlation functions needed 
in the analysis for ben ding moment. These quantities may be referred to 
as primitive cross-correlation functions corresponding to | r and | s . 

It is of utmost importance to keep the subscripts and superscripts in 
their proper order to avoid error. After substitution, the autocorrela- 
tion of bending moment reads 
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MgMsCx) 


N N M M 

im- 


i j r s 


r>“ r,s 

/ (-v)f ± fj(x - v) dv 

W —00 


Notice that the cross -correlation function fj_f j occurs in the expres- 
sion for bending-moment autocorrelation in the same form as in the accel- 
eration autocorrelation. Before proceeding to a discussion of this quan- 
tity in detail some remarks should be made about the analysis at its 
present stage of development. First of all, it is quite evident that 
the analysis -will reach vast proportions if the numbers N and M are 
very large. For example, if three rigid-body modes are used in conjunc- 
tion with three elastic modes then 36 primitive cross correlations are 
needed for the acceleration autocorrelations, while 324 primitive cross 
correlations are needed for the bending-moment autocorrelation. Some 
idea of how this goes may be seen in table I. 


Cross-Correlation Function f-jf j 

Consideration is now g iven to the determination of the cross- 
correlation function fj_f j keeping in mind that the analysis is becoming 

extremely complicated and that some simplifying assumptions may be made 
in the interest of retaining an analytical procedure which is -useful for 
practical applications. It is shown in the appendix that a typical term 
of the forcing-function matrix corresponding to the jth row can be written 


r b , . r 00 © 

£ - Xq - X 

/ Aj(y) d y / y 


J _b J 1 -00 

c (y) 


®i(x,y) dx + 




- Xp 1 - X 

c'(y) 


a* i (x,y) dx 


The function ijr may be identified as Kussner's function. The 
argument of the function as above written takes into account the effects 
of sweep and taper. The Kus sner function, as is well known, takes into 
account the unsteady aerodynamic effect of the lift buildup on the air- 
foil due to variations in <mj_. If the variations in a-j_ are slow compared 
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with the variation in ijr, over, say, 10 chords, then quasi-steady aero- 
dynamics may be used, and the quantity ijr may be replaced by a unit 
impulse. 5 in the following discussion it is assumed that Oj_ variations 
are slow enough to warrant this simplification. Under this simplifying 
assumption the forcing function reads 


pb p b * 

f jJ “ / d y Aj(y)<Xi(£ - Xo,y) + / dy A^'Cy)^^ - Xo*,y) 

-b -b 1 


Notice that t he e ffects of sweep and taper are still included. The cross 
correlation f 1 fj(X) may be computed in a straightforward calculation 

vising the above relation; namely, 


fjfj(X) = JJ dy dq Aj/yjAjCri) cy.(£ - Xq^ojXS + X - x^T]) + 

pb pb T ■ — ■ 

J J * dy drj A i (y)Aj , (n) cy.(£ - x 0 ,y)<Xi(£ + X - x 0 ',q) + 


p U pu 

/ / dy dq Ai'CyjAjCT}) ai(? - XoSyJcyXS + X - x^ii) + 

J _b J -b 1 

JJ ' dy drj Ai‘(y)Aj f (rj) OjU - x 0 , ,y)a i (C + X - x 0 ',n) 


If the turbulence is homogeneous, the quantity <Xj_a-j_ may be identified 

as a spatial autocorrelation function in two dimensions in which the 
spatial dependence is given in terms of the differences | y - 1} | , 

| xo - xo'j, and so forth. It Is important to notice that in the case of 

a swept wing (or tail) x Q is a function of the spanwise coordinate. 

In partic ular , if T is the sweep angle at. midchord then " 


xo(y) - xo(n) = sln r > (y - n) 


^Since a "frozen" turbulence field has been postulated, it may be 
stated equivalently that unsteady aerodynamic effects may be neglected 
in the cases where the turbulence integral scale is large compared with 
a characteristic dimension of the wing, in this case the semichord. An 
argument which leads to this simplification is given in reference 3. 
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Also, the last term on the right-hand side of the relation for fj_fj(x) 

is small compared -with the other terms, being approximately in the ratio 
of tail area to wing area squared. Therefore, this last term may be 
neglected with little error. Finally, if the autocorrelation for ctj_ 

is given the symbol cp^, the relation reads 


f i f j( x ) - JJ dq AiCyjAjCq)^ 

X + sin T(y - T]),|y - t] 

pb pb 1 

J J b . dy dn 

x + xo-xo‘,jy-n| 

pb pb 1 

J b J d y dT l A 1 '(y)A iJ (Ti)cp aa 

X + Xq' - Xq, |y - t] 


Notice that the effect of sweep is taken into account by including span- 
wise dependence in both coordinates of cp cax (x,y) . For the first integral 

this dependence is rather simply included, but in the last two integrals 
a more general relation must be 'used, namely, 


Xo(y) - x 0 '(ti) = l Q - sin I*,y + sin r,Tj 

where Z 0 is the distance between the midchord of wing and the midchord 

of the tail at the plane of symmetry. For most aircraft wherein both 
the wing and the tail employ moderate sweepback, the right-hand side of 
the above relation may be replaced by l Q . This aff ords considerable 
simplification in the last two integrals affecting f^f j. 

Certain special cases will now be considered in light of these 
results. The first case studied is for straight or moderately swept 
wings. Another special case will be considered which leads to somewhat 
simpler results as will be seen, that is, the case when cp txa (x,y) may 

be replaced by cp^Cx) . 


Special Cases 

Straight wings or wings with moderate sweep .- In order to avoid 
mathematical encumbrances, only o ne te rm will be considered on the right- 
hand side of the expression for f^f j remembering, however, that this 
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is only part of the total contribution. The remaining two terms should 
contribute in the order of 10 to 20 percent so that they should not be 
ignored completely. For this case, therefore, one may write 


w ■ £ dy dn A^yjAjCnteoaWy - ill) • • • 


Considering this expression in light of the definition of A (see 
appendix) , it is evident that an important special case is obt ain ed when 
i — j = 1 , that is, the case of the total lift on a rigid unswept (or 
moderately swept) wing restrained against motion. For this case, 


= qc(y)C la (y) 


Other special cases (e.g., rolling moment on a rigid restrained wing) 
may be obtained in a similar manner from A^. In order to evaluate the 
expression for f-jfj, a specification is needed with respect to tp^. 

In other words, the "space-time" autocorrelation (space-time correlation 
is directly proportional to space-space correlation under Taylor's 
hypothesis) for the turbulence is needed. Here, indeed, is the crux of 
the matter, but it is this quantity about which the least is known so 
far as atmospheric turbulence is concerned. The space-time correlation 
Is extremely difficult to measure even under ideal experimental test con- 
ditions. Furthermore, very limited information is available which affords 
a rational basis for estimating this quantity theoretically. 

If the turbulence is assumed to be Isotropic, then theoretical 
expressions could be used, but it would be difficult to justify their 
\ase. For example, an expression which is based upon isotropic turbulence 
considerations and which is not too difficult to handle mathematically Is 
the Gaussian curve 


^act.( x >y) = “(x 2 + y 2 )/L 2 

where w^ is mean- squared intensity and L is integral scale of turbu- 
lence. From the remarks it is clear that something more definitive is 
needed in order to complete a gust-load analysis of a full-scale aircraft, 
especially when spanwise partial correlation effects are to be included. 
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One -dimensional gust structure . - The case of the one-dimensional 
gust structure affords considerable simplification in two distinct areas 
of the problem, namely, (l) the determination of the primitive solution 
and ( 2 ) the determination of the autocorrelation. If the gust structure 
is one -dimensional the terms enteri n g the right-hand side of the dynamical 
equations of motion include the angle-of -attack f luctuat ions aj_ as a 

common linear factor. The matrix Fj j may be replaced by 



where the quantities f ^ and f j j are identical except that the angle- 

of -attack fluctuation does not erffcer into the spanwise integration. 

For this case, only one primitive solution is required for any particular 
quantity of interest. This is true, for example, in computing either 
center-of -gravity acceleration or bending moment or any other quantity 
which may be of interest. All that is required to arrive at a solution 
to a particular problem is to replace o^(£) by the unit impulse, to 

compute the desired primitive solution, and to use this primitive solu- 
tion in a superposition on 0^(5) to obtain the final desired result. 

Suppose, for example, it is desired to determine the center-of-gravity 
acceleration and the autocorrelation function of the wing bending moment. 
The primitive solutions may be designated and Mg*, respectively. 

Thus, 


Z(Z) 


Mg(Z) 



2*u - £)ai(£) d5 


m b *(z - S)a±(£) as 
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It is evident that this is a much simpler means of deter mining the final 
results. The analysis for this case is in all respects similar to the 
analysis of an electrical circuit with stochastic inputs. Accordingly, 
the relations which are developed in the analysis of such systems may 
be directly applied here. It is interesting to compare the above scheme 
with the previous one in terms of relative complexity. Such a comparison 
is best brought out by considering the number of primitive cross-correlation 
functions needed to complete the above analyses. In the case of the center- 
of -gravity acceleration it is necessary to compute only one pr imi tive cor- 
relation which is, of course, an autocorrelation rather than a cross cor- 
relation. In the case of the bending moment, it is similarly necessary 
to compute only one primitive autocorrelation. However, for bending moment, 
it is important to notice that the primitive solution depends upon all nor- 
mal modes | r *, and the total bending moment must be obtained from the 

relation 



Arl r * 


where the values of | r * are obtained by replacing 0^(5) by the unit 
impulse. 

The primitive -solution requirements and associated primitive- 
correlation- function requirements may be compared with those in table I. 

It is clear that the one-dimensional case is several orders of magnitude 
simpler to apply. One may generalize this result in the following way. 

For the one-dimensional case, let p represent some physical quantity 
of interest which depends upon the solution of the dynamical equations 
of motions. Let p* be the solution of the equations for cuj_( = 5 = 

Unit impulse. Then, the autocorrelation of p is given by 

<Ppp(x) “ f 'PpV*^ " £)<Paa(E) 

where is the autocorrelation of p* and q> nr t. is the autocorrela- 

tion of angle of attack, a. 
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If the power density spectrum of (i is of interest, then (ref. J, 

p. 322) 


$>(!(! (oi) = | M^ 'Wl : "®oa (ai) 


where $ is mean-power density spectrum and M*(cu) is obtained by letting 
a± = exp tcufc in the dynamical equations of motion. 

This relation is a familiar one which is basic to the study of sta- 
tionary time series and linear passive systems. It may be applied to 
any linear system with constant parameters provided the input quantity 
is a stationary random function of one independent variable. For a linear 
passive electric circuit these conditions are satisfied. For an airplane, 
these conditions will, in general, be satisfied only when the turbulence 
scale is "large." In this case the turbulent velocity fluctuations which 
influence the airplane response at any instant may be considered constant 
in the lateral (spanwise) direction. The analysis in these cases of large 
turbulence scale is a great deal simpler than the general analysis. Flight 
measurements of the power density spectrum of the gust considered as a 
one-dimensional process (refs. 3 > l6, and 17) may be used in a one- 
dimensional analysis of the type discussed herein. 

The function M*((u) appearing in the above relation may be cal- 
culated from the dynamical equations of motion. For example, if 70 points 
in the power density spectrum are to be calculated and if the analysis is 
for six degrees of freedom, it will be necessary to perform 70 inversions 
of sixth-order matrices to obtain M*(co) . 


SOME SIMPLE ILLUSTRATIVE EXAMPLES 


The general formulation of the problem of airplane gust response has 
been examined. This problem becomes exceed i ngly complex and difficult to 
handle analytically when instantaneous spanwise variations of the gust 
field occur. Several simplify ing assumptions were introduced into the 
formulation. But the assumption which achieves the greatest simplifica- 
tion in the analysis iB the postulation of a one -dimensional gust structure. 

For this case, the response may be characterized by the relation 


fcpp.Ctu) = |m/(co) [ ®aa(a>) 
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where 

p generalized dependent (response) variable 

$p.i(q)) power density spectrum of p 

Mu^ai) transfer function of p corresponding to sinusoidal gust 
input (obtained by solving dyn ami cal equations of motion 
for airplane) 

The quantity ,(<») is the power spectrum of w/u and may be 

represented by the one -dimensional isotropic spectrum (ref. J, p. 322 ) 
as 





where L is the integral scale of turbulence, and w 2 is the mean- 
square gust velocity of turbulent velocity fluctuations normal to the 
airplane flight direction and lying in the plane of symmetry. The mean- 
square lift response as a function of the parameters of the problem will 
now be investigated. In general, the transfer function M*(co) is a func 
tion of reduced frequency a£>cj U and certain nondimensional airplane 

parameters, for example, airplane relative density, pitch radius of gyra- 
tion, tail length, center-of -gravity position, tail effectiveness factor, 
and so forth. 

Using that inf ormation, one can write the following relation for the 
mean-square lift coefficient (ref. 18 ) : 



$ Cl C l (<x>) dn>= jT I M Cl *( 00) 2 ®oa(“) dco 


J | MCl*( £ b c/ L , K, l, It 


>•••)! 


yo as 
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where | = Thus, it is seen that the mean- square lift coefficient 

is a function of a number of nondimens ional airplane parameters and the 
parameter L/b c , which is the ratio of turbulence scale to wing semichord. 

The lift coefficient is also proportional to the turbulence mean-square 
intensity. Any other response variable of interest, for example, bending 
moment at wing root, would similarly depend upon these numerous parameters 
though not in the same way, of course. 

A complete parametric investigation (see, e.g., ref. 15) of the 
problem as above formulated is beyond the present scope of study. However 
it is desired to consider the following important questions: 

(1) What is the importance of unsteady aerodynamics in the lift 

buildup due to the gust? 

(2) What is the importance of the pitching degree of freedom? 

(3) How do the results compare with results obtained by the sharp- 

edged-gust formula? 


Importance of Unsteady Aerodynamics 

In a paper by Fung (ref. J) , the problem of unsteady aerodynamics 
was considered in the framework of a single-degree -of- freedom vertical 
translation. Results were obtained which included the effect of unsteady 
aerodynamics and also apparent mass effects but which neglected contri- 
butions of the horizontal tail. In the general formulation presented 
here the unsteady aerodynamic effects and apparent mass effects have been 
neglected, but the contribution of the horizontal tail has been included. 
It is desired to isolate the unsteady aerodynamic effect. To do this, 
consider the integral 


Cl 2 

( C L 2 )ref 


-r 

3t Jq 


| 2 (i + 3! 2 ) 


. n 2 


(SA +t a 

\b c kJ 

(1 + I 2 ) 2 


<U 


where k = 2m/itpSb c is the airplane relative density. This may be com- 
pared with the expression given in Fung's paper. 
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The reference lift coefficient is based upon the sharp -edged- gust 
formula 


(Cl) ref = Sjtv'/U 


where v 1 is the root-mean- square gust intensity. Figure 2 shows a 
direct comparison of results. It can be seen that unsteady aerodynamic 
effects become less important as the integral turbulence scale becomes 
large relative to the wing semichord. Also, for large values of airplane 
relative density, the effect of unsteady aerodynamics decreases for given 
values of L/b c . This latter influence is much less apparent than the 
scale effect. 


Importance of Pitching Degree of Freedom 

A parametric study of gust loads in two degrees of freedom was made 
by Brenner and Isaksen (ref. 18) . In that study, the gust structure was 
greatly simplified to the case of (l) a uniform step- function gust and 
(2) a uniform ramp-function gust. In what is perhaps the most exhaustive 
parametric treatment of the problem so far, charts and numerical data are 
given for more than 225 case studies. Certain effects on the gust response 
were investigated, for example, effects of airplane relative density, 
effects of pitch radius of gyration, effectB of tail length and tail 
effectiveness, and so forth. A parametric study of such scope is outside 
the present investigation. However, to settle the question of the Impor- 
tance of the pitching degree of freedom, a Similar analytical study for 
stochastic gusts is warranted. Such a program could make use of results 
obtained by Brenner and Isaksen if sufficient data on the asymptotic 
behavior of their responses are available. 

Only peak values of the responses were of interest in their study, 
so the responses for later times are not given in the report. To make 
use of step-function response data for computing responses to stochastic 
gusts, a Fourier integral operation is involved, and this integration 
extends over the range 0 <■ t ^ °o. Unfortunately, the data of refer- 
ence 18 do not cover the required range on t. 

For the sake of simplicity, two special cases are considered which 
axe characterized by two extreme values of the pitch radius of gyration, 
namely, J — »0 and J — >». The case for infinite radius of gyration is 
essentially the one-degree-of-freedom case of vertical translation, since 
with infinite moment of inertia about the pitch axis, the airpl an e can- 
not respond in pitch. For the other case, the pitch radius of gyration 
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(moment of inertia in pitch) is so small that the airplane offers effec- 
tively no inertial restraint about the pitch axis. Motion about the 
pitch axis is thus entirely the result of aerodynamic moments about the 
center of gravity. Any actual case would lie between limits given by 
these two simple cases. To establish these limits , the dynamic equations 
of motion must be solved corresponding to the two stated limits of pitch 
radius of gyration. 5his procedure leads to an expression for the 
transfer function in terms of the airplane parameters and the reduced 
frequency. If, in addition, one considers that the center of gravity of 
the airplane is subject to a narrow range of variation, say between 
25 and 50 percent of the chord, some additional, simplifications in the 
expression for the transfer function are possible. An expression which 
closely approximates the transfer function for small pitch radius of 
gyration is included in the following relations for mean-square lift 
coefficient with center of gravity at 50 percent of the wing chord: 



1 

jt 



|2(i + 5 |2) d| 



_ _ 

2 


2L 

+ | 2 


b c KZ’(l + R) 

< 


J 


(l + £ 2 ) 2 


Here, two new parameters have been introduced, the nondimensional tail 
length l ' , referred to wing semichord, and the ratio of tail area to wing 
area R. 


The ratio of C ^ 2 for J — >0 to Cj 2 for J— > » may be computed 
by referring to previous results, namely. 



r 

^0 


| 2 (1 + 3 £ 2 ) 


d£ 


[2L/b c wZ(l + R) 2 + | 2 (l + I 2 ) 2 


L 


£ 2 (l + 3£ 2 ) 


cL£ 


0 ( 2 L/boi«r + £ 2 (1 + i 2 )‘ 


This relation is plotted in figure 3 for the value Z(l + R) = 4.0. 
For the given parameters, the relative effect of the pitch degree of 
freedom becomes more important as the turbulence scale becomes large 
relative to the wing semi chord. Also, the effect of the airplane relative 



32 


NACA TN 3679 


density is to decrease the effect of the pitch degree of freedom as the 
relative density increases; all other parameters remain fixed. An increase 
in tail length will have the same relative effect as an increase in air- 
plane relative density. A larger tail will have a mitigating effect on the 
gust loads. 


Comparison of Results With Sharp-Edged-Gust Formula 

To see how these results compare with the sharp-edged-gust formula, 
some representative curves are shown in figure 4. These curves were cal- 
culated using the formula 


C L 2 _ 1 f 00 £ g (l + 3I 2 ) d| 

( Cl 2 ) ref * ^0 (2L/b c ic ) 2 + | 2 (l + | 2 ) 2 


These results show a decrease in mean-square lift response with 
increasing L/b c . Also, the effect of decreasing the airplane relative 
density is to mitigate the gust load (compare this with the relative 
effect of k on the loads when pitching is included) ; all other param- 
eters remain fixed. This latter effect is in accord with the result's of 
Brenner and Isaksen for step-function gusts. Perhaps .the effect of 
increasing the turbulence scale is tantamount to increasing the gust- 
gradient distance for ramp-type gusts . Bie overall effect on the lift 
response to stochastic gusts is in qualitative agreement with ramp-gust 
responses if one uses this reasoning. It must be remembered that this 
is a statistical process, and, if a root-mean-square load factor of, say, 
0.8 is predicted, this load is actually exceeded many times, as evidenced 
by the distribution of the load factor in the probability sense (ref. 19 ) . 
A load factor of 2.0, for example, is not impossible. In fact, If the 
distribution Is Gaussian, the probability of this happening Is 1.24 per- 
cent, and the odds against the load factor equaling or exceeding 2.0 are 
only about 80 to 1. Thus, while the sharp-edged-gust formula gives con- 
servative results for root-mean-square load predictions, it may be 
perfectly adequate for design. 


CONCLUDING RMAEKS 


An analytical scheme has been derived which can be used to estimate 
autocorrelation functions and, by use of a Fourier transformation, 
power density functions of airplane responses to atmospheric gusts. The 
analysis takes into account the principle degrees of freedom of the air- 
plane and includes effects of spanwise (as well as chordwise) variations 
in the instantaneous gust pattern. The final results depend upon the 
solution of dynamic equations of motion for the airplane for a set of 
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relatively simple forcing functions, namely, a set of delta functions. 

These solutions are superposed by convolution with, the act ual gust pat- 
tern to obtain the desired results . Two typical responses of interest 
are discussed, the center -of -gravity acceleration response and the wing- 
root bending -moment response. The autocorrelation functions of these two 
quantities are expressed in terms of a multiple summation, extending over 
the total number of degrees of freedom. Each term in the summation 
involves a correlation of the primitive solution cross correlation with 
a cross-correlation function. TOis cross-correlation function can be 
calculated by postulating that the turbulence field is stationary and 
homogeneous; also, the autocorrelation function in two dimensions must 
be known for the gust pattern. TOis quantity is defined by T^Cx ^) } 
it is by means of this function that the lateral gust field is brought 
into the analysis. Unfortunately, little is known about this function 
for the atmosphere, indicating that practical application of these results 
must await the determination of cp act (x,y) from atmospheric-turbulence 
research investigations. 

A special case which is vastly simpler in principle follows from 
neglecting the lateral variations in the gust pattern. This case has 
been labeled "one-dimensional 11 gust structure. The simplification of 
the analysis is such that the amount of work can be reduced by a factor 
of at least N^, if N iB the number of degrees of freedom. For this 
case, only one -dimensional autocorrelations enter; for example, cp aQ ,(x,y) 

may be replaced by 9 0UX (x) . Flight measurements of the power density 

spectrum of the gust considered as a one-dimensional process may be used 
in a one -dimensional analysis of the type discussed here. 

Some simple illustrative examples are treated for the case of a 
one-dimensional gust structure. The validity of neglecting the effect 
of unsteady aerodynamics is investigated by comparing calculations for two 
cases of one degree of freedom. The results show that this assumption may 
be justified only when the integral scale of turbulence is large compared 
with the wing semi chord. A corollary to this conclusion may be stated as 
follows: When conditions are such that unsteady aerodyn ami c effects are 

significant, then the lateral gust effects may become important also, and 
it does not seem consistent to include the one and neglect the other. Both 
should be either included or neglected in a rational analysis. 

In a one-dimensional analysis, it is conservative (i.e., higher 
loads are predicted) to neglect unsteady aerodynamic effects. It is 
not possible to state at this time the comparison for a two-dimensional 
analysis. 

The importance of the pitching degree of freedom is investigated 
by solving the two equations of motion (for heaving and pitching) under 
two extreme values of the pitch moment of inertia, namely, for a pitch 
radius of gyration of 0 and one of ». The latter case effectively 
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includes only the plunging degree of freedom. The results of this simpli- 
fied study show that the pitching becomes relatively more important as the 
turbulence integral scale becomes large relative to the wing semichord. 
Before one can draw general conclusions, however, an extensive parameter 
investigation must be made. 

The results of a single-degree-of -freedom analysis (neglecting unsteady 
aerodynamic effects) are compared with the results of the sharp-edged gust 
formula. For this comparison, the root-mean- square gust intensities are 
set equal. These results show that the root-mean- square loads are always 
less than those predicted by the sharp- edged- gust formula and become even 
smaller as the integral scale of turbulence increases. 

One final point is worthy of mention, namely, that this is a statis- 
tical process, and if a root-mean-square load factor of, say, 0.8 is pre- 
dicted, it should be remembered that this load is actually exceeded many 
times, as evidenced by the distribution of the load factor in the prob- 
ability sense. A load factor of 2.0, for example, is not impossible. 

In fact, if the distribution is Gaussian, the probability of this hap- 
pening is 1.2k percent, and the odds against the load factor equaling or 
exceeding 2.0 are only about 80 to 1. Thus, while the sharp-edged-gust 
formula gives conservative results for root-mean-square load predictions, 
it may be perfectly adequate for design. 

Extensive research investigations should be carried. out (l) to 
determine much needed information on the atmospheric gust structure, 
particularly on two-dimensional turbulence spectra so that a direct com- 
parison can be made with data now available and (2) to study methods of 
applying atmospheric turbulence data to practical problems of flight. 


Massachusetts Institute of Technology, 
Cambridge, Mass., July 1, 1955* 
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APPENDIX 

DERIVATION OF MATRICES ENTERING INTO DYNAMICAL 
EQUATIONS OF MOTION 


The equations of motion for an aircraft flying at mean velocity U 
through turbulent air are represented by the matrix equation 

[m] jq} +[a]jq| + [a]{ qj. + [e]{q} = {l} 


The various square matrices will now be considered in detail. 


Inertia Matrix Fm.1 


The inertia matrix may be derived by noting that the column 
matrix [m] jjgfj> follows from 


i r ? sv 

dt Jg Bq 


dm 


for q = z, 0, and |„(r = 1, 2, 5, . . . (in this order)) where 


V = z+x0+y$+^ A^Cx.y)^ 

r 

For example, one may consider the terms in the third row of [mjiql. 

These may be written ^ 

' at f s v | am ‘ f s ** ft [ 4 + + ^ + I 

= z / day + 9 / dmxy + 0 / dm y 2 + > V / dm A r (x,y)y 
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whence 

J ^ m 3 jj _ S x* ^xy> ^ 

where 


s x 

static moment at x-axis 

I xy 

product of inertia at xy-axis 

*x 

moment of inertia at x-axis 

J 

r\ 

dm A^x^yjy 
S 

11 

pj 

tP' 

dm A2(x,y)y 
a 


a 01' a 02 J ‘ * • 


and so forth. The determination of the remaining rows of [mj 
from an obvious extension of the above equations as 



= zM a + GSy + $S X + ^ ? r a zr 

= zSy + 01y + 01^ + V ? r a 0r 



{?} = 


= za^2_ + 9a, ( 


'91 


+ 


0°$! + £l M l 


where 




y 


mass of airplane 

static moment about y-axis 

moment of inertia about y-axis 


follows 
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Oar = jT dm A r (x,y) with r = 1, 2, 3, 


a 6r “ X dm A r (x,y)x with 


r = 2, 3 j • 


M T 


= J' dm A r 2 (x^y) 


It should, be noted, that use has been made of the orthogonality relation 
between normal modes in the above development. Also, depending on the 
use of either symmetrical or antisym me trical normal modes, certain terms 
will be zero. Furthermore, since the origin of the principle-axis 
coordinate system is here taken at the virtual center of gravity of the 
aircraft, all static moments S x , S y , and 1^ will be zero. For normal 

modes it may be shown that the terms oc^ r , a Qr , and a^, are zero also. 

However, they are retained here for generality. Thus, the mass matrix 
may be written 


m a 

% 

Sx 

“zl 

a z2 

% 

J y 

I xy 

a ei 

0*32 

Sx 

x 3cy 

lx 

a 0i 

“02 

“zl 

“ei 

°^1 

Mi 

0 

a z2 

°02 

“02 

0 

Mg 


which is symmetric about the main diagonal. 


Aerodynamic Matrix 



The aerodynamic matrix 



derives from those 


work expression due to 4* These can be traced to 
due to motion of the aircraft. A typical row from 



as obtained from dU e /B$ follows: 


terms in the virtual 

the aerodynamic forces 
this matrix, say 
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r 


o 

°3A 


w-y:: 


y <3y 


j -1,2,3,. 


^c(y)Ci (y) 

C8 


Z + ex + 0y + ^ ApCx^y)^ 


U 




^ + 


o --h 1 


y ay 


q.'c'(y)Cj '(y) 

Co 


z 4 - 9x7 + 0y + A^x^y)^ 


U' 




= z 


Z r\ u | r *^ 1 

- / ay yc(y)C z (y) + ^- / dy yc T (y)Cj ’(y) 
U J_ h a, U* J_ b r ‘a 


£j r»l3 . 1 r*^ ! _ 

- / dy yxc(y)Cj (y) + Jy / dy yx'c'(y)C l '(y) 

U o _>> ^ u v -"b ' 


0 


-r 


§ J dy y 2 c(y)Cj a (y) + §J- ^ ^ dy y 2 c'(y)Cj a '(y) 


z •» 1 / dy A r (x,y)yc(y)c z Jy) + 


7 ■ 

r / dy A r (x',y)yc’(y)Cj '(y) 

J _b' a 





+ 



r 
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where the definition of the elements a^, &J 2 ’ • • • is obvious. 

Continuation of this development yields the remaining matrix elements: 


j-1,2,3,... *- 


[hj] {?"} 3 J' h ** + b J b , ** c '(y)Ci 0 . , (y^ + J ^ ay =(r)c^Cy)5 + J ^ ay c'Cytej^'Cy)* 

1 § ^ + iP f _ b , 37 y='(y>c^'(y)j + X «r § /_ b 4 y 


:»y)c(y)C{ a (y) + 


b S , ^ ^^yWylo^Cy) 


jjaajJ ^4| - *[l / ®y “(yJO^Cy) + |r ^ ay £o’(y)Ci a ’(y) + e § ^ ay xxetyJCjJy) + §r ay £ , 5 r c'(y)Cj a '(y)l 

j-1,2,3,... L _t ■ J L 

§yy ay S&reCyJe^Cy) + ^ J ^ fly x'yc^yJCj^'Cy) + ^ j r | J ^ ay 2* P (x,y)c(y)aj a (y) + 

/ ay x'Ar(x',y)c'(y)Oj a '(y)J 


J-1,2,3,. 


{i}“ * § J ay A.i(x,y)c(y)Cj a (y) + ^ J ay A 1 (x',y)c'(y) c ^ x '(y) + B § yy ay xAxCxjylc'CyJOjJy) + 


pr / v ayx-A 1 (Sc,y)= , (y)c la -(y) + ji|§ ^ ay y*i(x,y)0j a (y) + y^ ay yA. 1 (x , ,y)e'(y)Ci a , (y) 


5 r | ^ fly AxCxjyjAjC^yJcCyJc^y) + §,• 3y A 1 (x',y)l r (x , ,y)c , (y)C la , (y) 


This matrix is not symmetric since x / x. 
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Aerodynamic Matrix jjx^] 

The aerodynamic matrix [a] may he derived in the same way as the 
matrix except that one considers only the static terms <j^q J, , which 

contribute to the "generalized force" cffi e /dq. For example,, the third 
row of [a] is given by the expression 


M{ 4 } = /. b - yay 

j— 2 , • • • 


qc(y)C z (y) 
a. 


L 


l 


0 + > e r (x,y)£. 




«b' 

r 

“ v __ 

•\ 

/ -y dy« 

/ 1 

q' c '(y) c z a '(y)j 

0 + 2^ ©r(xSy)l r 

y 


S- 

L — _] 

j 


= 0 


rib _ „b' 

q / -y dy c(y)Cj (y) + q 1 / -y dy c'(y)C7 ’ (y) 

J -h a J -b » a 


i •- * /: -y dy ® r (x,y)c(y)C la (y} + 


— r b ’ 

q' / -y dy 0 r (3T r ,y)c'(y)C z 1 (y) 

J _b’ a 


and so on.' The matrix JjxJ is not symmetric. 
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Elasticity Matrix 

The elasticity matrix [e] may be derived by noting that it derives 
from 8E^j8q_ where q = | and where 


E ± « 1/2 £ 


Thus, this matrix is a diagonal matrix of the form 



Forcing -Function Matrix J 

The forcing-function matrix does not depend in any way upon the 
generalized coordinates. However, as will be presently shown, it is a 
diagonal matrix whose elements depend upon the result of forming 8E e /dq 
for q = z, 0, . . . . 

Consideration may be given to that part of 8E e /ciq due to turbu len ce 
as follows: 

J ^ + 0y + ^ ^(x,y)| r qc(y)C 2a (y) J " TiyJojC^y) dr + 


a 

35 


J 1^2+®'+^+^ ArCX',y)£ r q.’c , (y)C^ a , (y) J ijr' (h - T,y)a i (r,y) dr 
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If the partials are formed, one obtains 


= 5 J c(y)C z ^(y) J *(t - T,y)a 1 (T,y) dT 


8 r ^ ~ r“ . 

g = i j <iy xc(y)Cj a (y) J ijr(t - T,y)a ± (T,y) dT + 


• • • 


d_ 

S0 


n b P 03 . 

= q J dy yc(y)Cj a (y) J t(t - ^yJo^C^y) dT + , 


•s _ ^b r°° 

- — = q. / dy A t ,(x,y)c(y)C 2 (y) / t(t - Tby^C^y) dT + . 

df J _■ b a ^ -co 


Using Taylor's hypothesis, it may be shown that a typical term 
from the above ensemble may be written 


5 " x o - x 




ct-i(x,y) dx 


where x = x_(y) = Distance parallel to x-axis from center of gravity 

O / 

to midchord and \|r = d^/dx. 

Note that Ut = £ and also that the turbulence ct i (x,y) is phased 
with respect to the aircraft center of gravity. Accordingly, the elements 
of Fjj are 


AA 


■s: 


dy A< 


r© 

\ - x o " x 

1 

8 

-€ 

0 

1 


(x,y) dx + 


p b* n“0 

/ <3y AS(y) / T|r 

J _b * J U -co 


S - x o’ - X 1 

c(y) 


Oi(x,y) dx 



HACA 03 ? 3879 


^3 

Ohe quantity Aj (y) must be chosen in accordance with the relations: 

A x (y) = ic(y)C z (y) 
b a 

Ag(y) = qxc(y)Cr (y) 

Ur 

A^(y) = qyc(y)Cj a (y) 

Ai,.(y) = q.A 1 (X,y)c(y)C^(y) 

A 1 '(y) = q'e’CyJC^' (y) 

Ag' (y) = q'xc’CyJC, ' (y) 

c ‘’a, 
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TABLE I.- DETERMINATION OF NUMBER OF PRIMITIVE CROSS 
CORRELATIONS NEEDED IN ANALYSIS 


Number of primitive cross 

correlations for - 

Center-of -gravity 
acceleration 

Wing -root bending moment 

N 

M 

N 

M = 1 

M = 2 

M = 3 

M = 4 

2 

4 

2 

4 

16 

36 

64 

3 

9 

3 

9 

36 

81 

144 

4 

16 

4 

16 

64 

144 

256 

5 

25 

5 

25 

100 

255 ■ 

400 

6 

36 

6 

36 

144 

324 

576 
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Figure 2.- Comparison of mean-square lift response using unsteady and 

quasi-steady aerodynamic theory. 







Figure 3*- Effect of pitching, degree of freedom on mean-square lift 
response. I ' (1 + R) = 4.0j center of gravity at 50 percent chord. 









